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ABSTRACT 


Pseudospectral methods are used to effieiently solve complex, non-linear optimal control 
problems and can be used to develop rapid and robust control solutions in astrodynamics 
applications. With the renewed efforts in resuming future manned missions to the moon 
or near-earth asteroids, fuel optimal flight trajectories are desired in order to alleviate the 
need for overestimates in fuel carried onboard, which has a negative impact on the 
overall spacecraft design. 

This thesis uses a pseudospectral optimal control solver, DIDO©, to solve for a 
fuel optimal moon-to-earth trajectory using only the auxiliary engines of the Orion Crew 
Exploration Vehicle (CEV) for the entire mission. Additionally, higher-order necessary 
conditions are examined to test the extremality of the computed solution. Singular arcs 
that appear in a main engine variable thrust, fuel optimal trajectory are examined by 
applying the Bellman Pseudospectral method. The feasibility of using the CEV auxiliary 
engines in place of the main engines for the singular arc maneuver is also investigated. 
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I. 


INTRODUCTION 


A. BACKGROUND 

The NASA Orion Crew Exploration Vehiele (CEV) program has more stringent 
safety requirements as compared to the Apollo missions of the 1960s and 1970s. One 
such requirement is for “anytime return” of the spacecraft in the event of a mission abort 
scenario. The increased safety margins that need to be incorporated into the spacecraft 
design begin to form a competing set of requirements. The design of spacecraft, whether 
manned or unmanned, is primarily driven by the mass of fuel carried onboard to perform 
maneuvers while performing its mission. This, in turn, drives the total mass of the 
spacecraft, which further constrains the design for mission essential equipment and can 
heavily influence the design and/or selection of the launch vehicle. Excess safety 
margins should be avoided, however, in order to maximize mission effectiveness. 

The application of optimal control theory to astrodynamics applications allows for 
designers to generate engineering solutions that maximize desirable performance 
characteristics as well as minimizing the need for excess design margins. Current design 
practices may unwittingly limit the solutions being created by enforcing seemingly 
practical constraints. As the spacecraft design matures, clearly there is a need to start 
imposing design constraints due to the available technology, product specifications, and 
material. However, early in the design it is possible to find “best” solutions by limiting 
the number of constraints imposed on the design and mathematically validating the 
solutions using existing optimal control theory. The tools for solving complex optimal 
control problems allow for rapid and robust solutions and straightforward methods are 
available for ascertaining their feasibility. 

B. SCOPE OF THE RESEARCH 

This thesis investigated the feasibility of using solely the lower thrust auxiliary 
engines of the Orion CEV for return to earth from low lunar orbit. Using DIDO©, which 
utilizes MATLAB to solve complex, non-linear optimal control problems, solutions were 
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generated whieh demonstrate the existence of trajectories which would allow the 
spacecraft to return using only the auxiliary engines while meeting specific mission 
timeline and fuel consumption requirements. Furthermore, the solutions generated are 
shown to meet the necessary conditions of optimality by applying Pontryagin’s Theory. 
While there exists a range of possible solutions, two significant solutions are presented; 

(1) the absolute minimum time required for the maneuver with no constraint on fuel, and 

(2) the minimum fuel required to reach the given final boundary conditions within 48 
hours. 

The first case is shown primarily to show the absolute lower limit on time 
required for the maneuver with the given engine characteristics and the dynamics of the 
problem. The necessary conditions for optimality of this solution are shown and a 
feasibility analysis is conducted to check the solutions. The feasibility of the control 
solution is conducted by the application of the Bellman technique. While the Bellman 
technique is not shown for the second, minimum fuel case, it demonstrates a method to 
verify the fidelity of the generated solutions and determine if the solution is “fiyable” 
given a nominal set of discrete control points. The second case shows that the auxiliary 
engines can be used to complete the entire moon-to-earth Trans-earth Injection (TEI) 
maneuver sequence and still meet the 48-hour time requirement and use less fuel that the 
capacity of the Orion’s fuel tank. In all cases, the necessary conditions for optimality are 
demonstrated using Pontryagin’s Hamiltonian Minimization Condition (HMC) and where 
singular arcs appear in the trajectory, further analysis is performed to show optimality. 

Another problem investigated in this thesis is the control of the singular arc 
segment of a previously determined DIDO© solution using the CEV main engines. The 
main engine solution resulted in the return trajectory consisting of three separate TEI 
bums. The first and third burns are maximum thmst burns, known as bang-bang 
maneuvers. The first burn raises the apoapsis of the low lunar orbit, while the final bum 
provides the necessary velocity to escape the moon’s gravity on its way towards earth. 
The middle bum is given as a finite burn of varying thmst where the spacecraft conducts 
a plane change maneuver. A study of the high thmst singular arc solution was conducted 
in order to attain more in-depth details of the control trajectory in order to determine the 
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feasibility of using the auxiliary engines to conduct the singular burn segment of the high 
thrust return trajectory. Figure 1 shows the thrust profile for the minimum fuel, fixed 
time return solution using the Orion CEV main engines found by Yan et al. (Yan, Gong, 
Park, Ross, & Souza, 2010). 
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II. EQUATIONS OF MOTION 


A. INTRODUCTION 

For interplanetary trajectories, the two body problem no longer suffices as the 
spacecraft will pass through areas of varying gravitational influence as it gets further 
away from one body and gets closer to another. Even outside the sphere of influence of a 
particular body, the gravitational acceleration of that body may provide a minimal, but 
influential gravitational affect, often referred to as a gravitational perturbation. This 
chapter first defines the specific reference frame used for each case in this thesis. 
Secondly, it describes the gathering of ephemerides data on the gravitational bodies 
which have significant effect on the spacecraft trajectory which is ultimately used in 
setting up the dynamics of the moon-earth return mission. Finally, a brief description of 
the dynamics is given, followed by an example of why the problems require n=4 for 
moon-earth return missions. 

B. DEFINING THE REFERENCE FRAME (J2000) 

If one has to move from one specific location to another, one has to have good 
knowledge of the current position as well as an accurate description or knowledge of the 
intended final position. In common terms, it is called navigation. On the surface of the 
earth, the art of navigation is generally simple. Barring any obstacles, one simply takes 
the shortest distance between two points (which are fixed) and travels along the line that 
connects the two points. This is done by establishing a reference frame. Assuming that 
the travel is in only two dimensions, a planar reference frame can be used. Longer 
distances on the surface of the earth can be referenced to the surface of a sphere. Despite 
the rotations of the earth, the two dimensional traveler need only stay on the line 
connecting the two points. The reference frame, whether a flat plane or surface of the 
sphere is considered inertial since it is not accelerating with respect to the traveler. 

Depending on the specific application, such as interplanetary space missions or 
tracking satellites in orbit, a number of different reference frames are used. For rough 
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calculations, the reference systems used ean use an objeet’s eenter of mass as the eenter. 
However, for spaeeflight in the vieinity of the earth and moon, it might be useful to use 
the earth-moon barycenter, or system eenter of mass. Synodie systems are also used 
whieh are systems that rotate about the baryeenter. In synodie eoordinate systems, it is 
assumed that the rotation about the baryeenter is with eonstant veloeity with respeet to an 
inertial frame. 

The Heliocentric Coordinate System (XYZ) is a sun eentered system with the 
primary axis, X, pointed at the vernal equinox. It is a right handed system with respeet to 
the direetion of earth’s orbit around the sun with Z perpendicular to the eeliptic plane. 
With the correet transformations, this eoordinate system could be converted into a 
synodic eoordinate system, by using the barycenter of the solar system whieh is a point 
slightly off from the geometrie center of the sun. 



Figure 2. Helioeentrie Referenee Frame (From Vallado, 1997) 


A eommonly used eoordinate system in astrodynamies is ealled the Geocentric 
Equatorial Coordinate System (UK) which is a non-rotating system with the primary 
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axis, /, pointed at the vernal equinox. The J axis is 90 degrees to the east and the K axis 
points out of the earth’s north pole. This system is also known as the Earth Centered 
Inertial (ECI) system. The ECI system in aetuality does move over time due to effeets 
sueh as precession and nutation. The J2000 system is essentially the ECI system at a 
fixed epoeh, previously noted. This is eonsidered an inertial frame from whieh 
caleulations known as reduction formulas can be used to advance (or regress) to another 
epoch. 



Eigure 3. Earth Centered Inertial Reference Erame (Erom Vallado, 1997) 

The problem with space travel is in determining and characterizing an inertial 
reference frame, since by nature all the objects in the solar system are in a state of 
continual motion, and in most cases non-uniform motion due to the effects of 
gravitational forces between them. Eor short time spans, it may suffice to assume an 
inertial frame with the earth at its center; however the nonuniform motion of the earth 
and other solar system bodies have an increasingly perturbing effect over time. 
Therefore, to accurately determine the position of a body (celestial body or spacecraft) 
traveling through the solar system, one has to have a method to accurately determine the 
locations of all the gravitational bodies, or at least the ones which have the most 
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significant effect on the motion of that body. The positions of the solar system bodies 
ean be determined by observation and the future positions ean be determined, assuming 
the motions are well eharaeterized, at a future time. 

For eonsisteney, it has been useful to designate a speeifie time, or epoeh, as a 
referenee time with whieh the referenee frame used for navigation purposes ean be 
advaneed or rotated. The eurrent system in use today is ealled J2000. The J2000 epoeh 
is speeified to be, January 1, 2000 at 11:59:27.816 TAI (International Atomie Time), or 
about noon UT (Universal Time). Therefore, J2000 refers to a speeifie epoeh from whieh 
to attaeh a eoordinate system as a basis for all eelestial observations (Vallado, 1997). 
Typieally, an ECI referenee frame is used as the basis. If a moon-eentered inertial frame 
was desired, then with the proper transformation, the frame eould be translated from the 
earth to the moon’s loeation at the J2000 epoeh. 

In addition to planetary motion about the sun, the Earth’s position and orientation 
with respeet to the J2000 ECI frame varies due to gravitational interaetions with other 
bodies in the solar system, but primarily due to the Moon and Sun. The four primary 
transformations whieh must be made are for precession, nutation, sidereal time, and 
polar motion. Preeession results from perturbations from the Sun, Moon, and planets and 
it ehanges the orientation of the eeliptie with respeet to the J2000 frame. The reduetion 
formulas for preeession allow for the ealeulation of the mean equator of date from the 
mean equator at epoch (J2000). Nutation is largely eaused by the moon and is the 
lengthiest transformation, eonsisting of over 100 trigonometrie terms and is a periodie 
perturbation. Taking into aeeount the effeets of nutation transforms the mean equator of 
date into the true equator of date. The third transformation, sidereal time, transforms the 
non-rotating true of date frame to the Earth-fixed eoordinate system. Einally, polar 
motion aeeounts for the ehanging loeation of the North Pole. The motion of the North 
Pole follows a eireular spiral pattern and has a maximum variation of about only 9 meters 
in any direetion (Vallado, 1997). 
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Figure 4. Precession and Nutation (From Vallado, 1997) 


The problems studied in this thesis use a J2000 moon-centered frame of reference 
since the spacecraft originates in close proximity to the moon. In this way, the moon is 
the primary gravitational body exerting a force on the spacecraft, while other celestial 
bodies, including the earth, exert a perturbing force on the spacecraft. 

C. JPL HORIZONS DATA 

The Jet Propulsion Lab (JPL) HORIZONS on-line system provides accurate 
ephemeris data for solar system objects to include 538186 asteroids, 3066 comets, 170 
planetary satellites, 8 planets, the Sun, LI, L2, select spacecraft, and system barycenters 
(Jet Propulsion Laboratory). Using the web interface, ephemeris data was collected for 
earth, sun. Mars, and Venus as target bodies with the origin body as the moon. The 
ephemeris data was collected in the form of position (X, Y, Z) and velocity (VX, VY, 
VZ) vectors from the moon to the target bodies in the J2000 frame covering a time span 
of ten days starting at midnight on April 2, 2024, until midnight on April 12, 2024, in 
increments of one minute. The vectors were with respect to the ecliptic and mean 
equinox of the J2000 epoch. 
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In order to have access to accurate position vectors to the target bodies at any 
time, the relative motion of the target bodies can be characterized as a smooth function of 
time. In this study the Matlab function, polyfit, was used to be able to compare the 
accuracy of the data used by Yan et ah, who used the same function (Yan et ah, 2010). 
Even though the data points are very close together in comparison to the total duration, it 
is desired to create a smooth function in order to determine the position vectors between 
the discrete points. The goal for the overall project is to ensure that the accuracy of the 
position of the spacecraft and the celestial bodies is within one kilometer. At spacecraft 
interplanetary velocities, one minute is actually a long duration. Take for example a 
spacecraft travelling at a velocity of 7.7 km/sec. In one minute, the spacecraft has 
travelled 462 km. 

The polyfit function in Matlab interpolates between the data points by fitting a 
polynomial of specified order between the data points, x, by solving for P(x) in a least 
squares approximation. In order to minimize calculation errors when solving for the state 
vectors using polynomials, the order of the polynomial should be minimized but balanced 
against desired accuracy. For the position vectors, an accuracy on the order of one meter 
was desired and was accomplished by using polynomials of degree n = 12. The 
coefficients for the position and velocity vector polynomials are given in TABLES 1 and 
2. The position and velocities of the bodies can then be determined simply by applying 
the following equation, where P applies to either the Earth or Sun position or velocity at 
the specific time. 

n 

P = ^ ,t = days from start of epoch (1 • 1) 

k=0 

Figures 5 and 6 show the polynomial fitting errors for the sun and earth positions 
taken from the Horizons ephemeris database; increasing the order of the polynomial did 
not result in an increase in accuracy. 
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Table 1. Moon-to-Sun Coefficients 
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Table 2. Moon-to-Earth Coefficients 
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Figure 5. Sun Position Error wrt Moon 



Julian Day (00^=2455347.5) 

Figure 6. Earth Position Error wrt Moon 
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D. THE IMPORTANCE OF ACCURACY 


One might ask why it is so important to take into account what seems like very 
small variations and perturbations. Obviously, over long periods of time, the effects of 
the perturbations would increase, but these could be taken into account using frequent 
observations and updates to spacecraft navigation systems and simply make course 
adjustments during flight. The greatest constraint on spacecraft is typically its mass, due 
to launch constraints imposed by launch vehicles. A portion of this mass is fuel which 
the spacecraft requires for orbit transfers, station keeping, and attitude control. While the 
spacecraft could conceivably make correction maneuvers in flight, it comes at a cost of 
fuel. Any additional fuel to a spacecraft comes at a cost to allocated mass in other 
systems. For interplanetary missions, especially for future manned exploration missions, 
the payload and support systems are critical to maximize mission duration and ultimately 
mission viability. In practice, the spacecraft design should minimize the amount of fuel 
required to conduct its mission. Therefore, for design purposes, it is important to take 
into account additional gravitational effects that act on the spacecraft. 

If the target body’s motion, such as the Moon or Mars, is very accurately 
determined, then optimal trajectories and thrust regimes could be used without the need 
for correction maneuvers during transit. An analogy might be a bullet travelling 
downrange to a target at a distance of 1,000 meters. If the bullet’s initial trajectory is off 
by 0.1 degrees, it will only miss the center of the target by 1 millimeter. However, 
suppose a spacecraft travelling in a straight line, travels a distance of 80,000,000 
kilometers, the rough straight line distance between Earth’s orbit around the Sun and 
Mars’ orbit. A trajectory error of the same, 0.1 degrees will have the spacecraft miss its 
intended target by over 120 kilometers. Considering the great distance travelled, this 
may not seem like a significant number, however when trying to hit the target on the 
mark, with minimum fuel, this is indeed a considerable error. This highlights the need 
for a very accurate description of the motion of the planetary bodies and this primarily 
relies on an accurate accounting of time from the inertial reference frame time. 
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For a spacecraft returning to Earth from the Moon, in order to minimize fuel, 
optimal control trajectories are being created such that the spacecraft has a 1 km window 
for earth orbit insertion. This means that in order for the trajectory to be accurately 
designed, the spacecraft dynamics need to be of very high order of accuracy. It can be 
shown that for a moon-earth return trajectory, the perturbation effects of the Earth and 
Sun must be taken into account in order to achieve the accuracy required. In order to 
demonstrate this fact, a portion of the trajectory from Yan et al. moon to earth trajectory 
using the CEV main engines was investigated to show the effects of changing the number 
of bodies included in the dynamics equations (Yan et ah, 2010). The segment 
investigated has the initial position immediately following the first TEI maneuver and the 
target position immediately prior to the second TEI maneuver. The trajectory was 
propagated using the MatEab ode45 function and allowed to propagate until the time of 
the second burn maneuver. Starting with a simple restricted two-body dynamics equation 
using only the Moon’s gravity (Body 2), successive runs were conducted, each adding the 
gravitational effects of the Earth (Body 3), Sun (Body 4), and finally Mars (Body 5). The 
time-varying positions of the earth and sun with respect to the moon were taken into 
account as described previously as well as the governing equations of motion as 
described in the following section. The initial conditions are as follows: 


^0 


406.4091528635830 

km 

3^0 


255.0671214283607 

km 

^0 


-1,824.834644517191 

km 



1.445606181317849 

km/sec 



-1.651932481918282 

km/sec 

3oJ 


-0.274681204800144 

km/sec 


Eigures 7 and 8 show that the four-body dynamics, which includes the moon, 
earth and sun, satisfy the accuracy requirements while the addition of the fifth body do 
not contribute a significant effect on the trajectory, therefore its effect need not be 
included. 
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Figure 7. Effects of Additional Gravitational Bodies on Displacement 



Figure 8. Effects of Additional Gravitational Bodies on Velocity 
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E. EQUATIONS OF MOTION 


I. The Restricted Two-Body Problem 

In the restricted two-body dynamic model, the mass of the spacecraft is deemed 
negligible and equation of motion is: 


r 


GM (f'^ 


\r J 


( 1 . 2 ) 


This equation describes the force of attraction between two bodies, namely a 
planet (moon or sun) and a spacecraft where G is the universal gravitational constant, M 
is the mass of the large body, m is the mass of the spacecraft, and r is the position vector 
indicated by the negative sign to be from the spacecraft to the large body (Figure 9). For 
generality, the origin of the inertial reference frame will not be located at the center of the 
large mass. Therefore, the position vector r will be the difference between the vectors 
from the origin to each of the bodies as follows: 


r = r — r 

' 'sc 'm 


(1.3) 
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m 



Figure 9. 2-body Dynamic Model 

Since the data from HORIZONS is in XYZ vector format, the r vector is 
converted into Cartesian coordinates, substituting /u^ioxGM . Three equations of 
motion are generated, one for each axis of motion where; 


x = - 




{r + y^ + z^),lx^ + r + t 

Muy 


z = - 


[x^ + y- + e)4x- + y^ye 

_ IV _ 

(x^+y-+z-)4x-+y^ + z^ 


- i, X x^^ x^ 


j, y = ysc-yM 


(1.4) 


^ ^SC 
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2 . 


The n-Body Problem 


From the point of view of an inertial observer, the restrieted n-body equation of 
motion with respeet to the primary body is shown in Equation (2.5) below. Figure 10 
shows a diagram of the n-body geometry when there are only three bodies (n=3). 


/“if 


1,2 


1,2 


■Z/', 

/=3 


^ f 


/,2 


V /.2 


+ ■ 


J 


(1.5) 
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Third Body 


ri,3 
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/ 
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/ 

/ r3,2 
/ 

/ 


■=^m2 


Spacecraft 


Figure 10. Three-Body Geometry (After Vallado, 1997) 

The first term in Equation (1.5) is the 2-body part of the equation of motion, 
where the subseript “1” represents the primary gravitational body. The spaeeeraft in the 
equation is body number “2.” The summation term on the right represents the 
gravitational perturbation forees due to additional bodies j through n, starting with body 
number “3.” The perturbation term itself has two parts, the left term ealled the direct 
effect beeause it represents the foree aeting direetly on the spaeeeraft by the body. The 
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right hand perturbation term is ealled the indirect effect beeause it represents the 
gravitational effeet on the primary, or number one, body. This form of the equation of 
motion for the n-body problem is known as the relative form and is used beeause the 
motion of the spaeeeraft is ealeulated relative to the primary body (Vallado, 1997). 

The indireet effeet of the perturbation effeets ean be negleeted if the motion of the 
n-bodies is adequately eharaeterized. Beeause the relative positions between the 
gravitational bodies ean be found by eelestial observation, the gravitational forees aeting 
between the bodies ean be ignored. Using the available JPL Horizons data, it is assumed 
that the effeets of the gravitational bodies on eaeh other are already taken into aeeount, 
therefore the only perturbation effeets that need be eonsidered aeting on the spaeeeraft 
are the direet effeets. Therefore, Equation (1.5) reduees to: 


r, ^ 




1,2 


i=3 




( 1 . 6 ) 


3. Equations of Motion for Moon to Earth Return Mission 

The speeifie equations of motion used in this thesis are driven by the 
establishment of a J2000 moon-eentered, translating Cartesian frame. The origin is fixed 
to the moon’s eenter and the primary axes are aligned with the J2000 sun-eentered 
inertial Cartesian frame (SXYZ). The SXYZ frame has the X-Y plane aligned with the 
plane of Earth’s orbit and the X direetion points to the Eirst Point of Ares. The SXYZ is 
a right handed system, with the Z axis perpendieular to the X-Y plane. 

Cartesian eoordinates are seleeted first, for ease of using JPE Horizons ephemeris 
data, and seeondly beeause it avoids singularities in spherieal eoordinates that may arise 
due to trigonometrie eonstraints (Yan et ah, 2010). 
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Figure 11. J2000 Sun and Moon Centered Frame (From Yan et al., 2010) 

The position of the spacecraft is represented by , where is 

the position of the Moon with respect to the center of the Sun, and fp^^ is the position of 

the spacecraft with respect to the Moon. Taking into account that the position of the 
spacecraft is with respect to the moon in a sun-centered inertial frame, then the following 
equation holds for fp /^; 


(1.7) 



Defining the velocity as; 


( 1 . 8 ) 



Taking the first and second derivative of Equation (1.7) with respect to the J2000 
Sun-centered inertial frame yields; 
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(1.9) 


rpis = ^M/s + + V/, + X rp,M 


'PIS 


^MiS "*" "*" "*" ^z, ~*~ ^MIS ^ ^PIM "*" ®M/5 ^ fpi^^ 


( 1 . 10 ) 


Noting that = 0 because the orientation of the Moon-centered frame is 
aligned with the Sun-centered frame, Equation (1.10) simplifies to: 






( 1 . 11 ) 


The fully expressed spacecraft dynamics also have to be expressed in the Sun 
centered inertial frame and will include the dynamics as in Equation (1.11) as well as the 
controls or forces that the spacecraft will generate as part of the complete equations of 
motion. There are three control variables which are necessary for describing the three 
dimension Cartesian components of the thrust. T is the magnitude of the thrust and a and 
P respectively define the azimuth and elevation angles. 

The dynamics can be succinctly expressed by applying Newton’s second law in 
vector form as follows: 


fnrp,s=G^+Gp+G,+T ( 1 . 12 ) 

Where m is the mass of the spacecraft, f is the thrust vector, and G^,Gp, 
andG^ are the gravitational forces of the Moon, Earth, and Sun respectively. The Moon’s 
translating acceleration, also a required portion of the overall system dynamics can be 
defined as: 
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(1.13) 



The thrust vector is defined in the Moon centered frame defining a as the azimuth, 
and P as the elevation by the following; 

T cos « cos fi 

TsmacosP (1-14) 

T sinfi 

The time-varying components of the displacement between the spacecraft and the 
center of the Moon are given by (x, y, z). Between the perturbing bodies of the sun and 

earth, and the primary body (moon), they are (X 5 ,y 5 ,Z 5 )and(x£,y£,Z£)respectively so 

that the magnitude of the displacements of the spacecraft from the gravitational bodies 
can be defined by: 




With representing the gravitational constants for the Moon, Earth, 

and Sun respectively and representing the exhaust velocity of the spacecraft, the final 

form of the system dynamics can be written, which includes the time-varying mass of the 
rocket in the following equations of motion. Equation (2.16) becomes the governing set 
of equations of motion that are used in the problems solved in this thesis. 
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III. OPTIMAL CONTROL THEORY 


A. THE OPTIMAL CONTROL PROBLEM 

Problems concerning “optimization” have a long history dating back to the 
ancient Greeks where through geometry they were able to determine answers to problems 
such as the shortest distance between two points (a line) and the maximum area enclosed 
by a given length of perimeter (a circle). These two problems are part of a class of 
problems known as optimization theory, where the former is a minimization problem and 
the latter a maximization problem. The calculus of variations is a theory that emerged 
during the eighteenth century that deals with the optimization of integrals. John 
Bernoulli (1667-1748) posed a problem concerning the minimization of an integral in his 
famous brachistochrome problem (Greek, brachist = shortest, chromos = time) of 1696 
involving a bead sliding under gravity along a smooth wire in which he asks what shape 
the wire must be to reach the end point in minimum time. During the 1950s, the theory 
of optimal control emerged as a direct consequence of space exploration efforts by the 
Americans and the Soviets. The theory would allow a spacecraft’s trajectory to be 
controlled in such a manner as to reach a desired destination using minimal fuel and in 
minimum time (Pinch, 1993). 

Optimal control theory, however, is not limited to the physical guidance and 
control of a spacecraft. It has applications for mission planning and in spacecraft design. 
Consider the case for a manned spacecraft mission in which time is a critical design 
factor. The mission duration drives requirements such as the fuel required for navigation 
and control of the spacecraft and amount of life support equipment and supplies needed 
to sustain the crew. Once the mission objectives are specified and the physics, or 
dynamics, of the problem are understood, the designer will seek for feasible solutions to 
problems such as the minimum time to be able to complete the mission and how much 
fuel to carry can be solved. Therefore, optimal control theory can be used in the design 
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process to determine feasible solutions to engineering problems. The strength of the 
results gained lie in the faet that the theory has a sound basis in mathematies, and ean 
therefore be proved in a straightforward manner. 

B. DERIVATION OF THE NECESSARY CONDITIONS USING THE 
CALCULUS OF VARIATIONS 

The problem being solved is to determine the eontrolw(t), sueh that the path, or 
trajeetory, from an initial state x(to) to a final state minimizes the some eost 

funetional / [x], where; 


7[x] = j ^ F(x,M)<it + £'(x^ ) 


( 2 . 1 ) 


The integral term represents the running eost and E^x^^) represents the end-point 

eost. Sinee real systems are being eonsidered, it is safe to assume that any trajeetory 
ofx(t) will be eontinuous and twiee differentiable within the domain of the state spaee. 
It is also allowable to have the eontrols be pieeewise eontinuous and bounded within 
some finite eontrol spaee. In order to show that a ourvex = x*(t) is optimal, neeessary 
eonditions ean be shown, under some assumptions of smoothness, using the ealeulus of 
variations and are as follows (Pineh, 1993): 

Hamiltonian 

H(A,x,u) = F(x,u) + f (x,u) 


Euler-Lagrange Equation 



du 


(2.3) 
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Adjoint Equation 



(2.4) 


Transversality Condition 

(2.5) 

Where, E{v^,x[tj^:= E{^x[tj^ + vJis ealled the Endpoint Lagrangian. 


The variable, A , represents the eo-state and while it has no physieal meaning, it 
has eorresponding values that mirror the state variables. The adjoint equation represents 
the eo-state dynamies and has a direet relationship to its eorresponding state variable via 
the Hamiltonian equation. The transversality eondition looks at the endpoint eost based 
on the final state in order to find final eonditions of the eo-states. 

C. PONTRYAGIN’S PRINCIPLE 

In 1955, it was diseovered that a problem existed with the Euler-Langrange 
equations in that it was impossible to solve for bounded eontrols. The problem was not 
with the engineering or the physies of the eontrols in question, rather with the math. 
Pontryagin then invented a new theory known as the Hamiltonian Minimization 
Condition (HMC) whieh then replaeed the Euler-Eagrange equations. The HMC was able 
to handle bounded eontrols, whieh is typieally the ease in any eontrol system. In most 
oases, steering or attitude eontrol is unbounded sinoe all the angles in a oirole oan be 
desoribed in positive or negative multiples of In . In the ease of an applied torque or 
thrust for eontrol, the realities of physios diotate that there will be upper and lower limits 
to suoh a eontrol. Therefore, Pontryagin’s HMC was therefore a neoessary development 
in the praotioe of optimal eontrol theory (Ross, 2009). 

An example of the oonstrained eontrol problem oan be desoribed in a two- 
dimensional ease of an orbit transfer problem where the state and eontrols are given by 
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x = [x,y,v^,v^,mf and u = [T,/3f. The state is deseribed by the position and veloeity 

in two dimensions, and mass. The eontrol variables deseribe the magnitude of thrust and 
angle with respeet to the referenee frame. The Hamiltonian for this problem beeomes: 


H{A,x,u) = A^v^ + Av K^{t) + —cos/3 +A K (t) + — sm/3 

V m J \ m 


+ / 1 ,. 




V 


( 2 . 6 ) 


The portion of the dynamies equations that do not depend on the eontrols are given 
by^, (t). The exhaust veloeity, v^, is eonstant and is a funetion of the physieal 

eharaeteristies of the roeket motor. The dynamies portions of the equation are exeluded 
for brevity beeause they do not eontain any eontrol eomponents. 

From the Euler-Lagrange equation, the optimal steering law is determined by: 
dH ^ , KyT ^ ^ 

-= —^smpH—^eosp = 0 (2.7) 

dp m m 

tany0 = i (2.8) 

^vx 

However, determining the optimal thrust law proves problematie sinee T vanishes: 

■^^ = —eosy0 +—siny^- —= 0 (2.9) 

dT m m 


The Euler-Lagrange equations give areas where the slope of the Hamiltonian is 
zero, thus giving potential loeations where maxima and minima may exist. However, if 
the slope is not zero at the endpoints, then these points would be exeluded even though 
they might also eontain maxima and minima. Example Hamiltonian eurves are given in 
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Figure 12 where the Euler-Lagrange equations would give an ineorreet value of the 
eontrols if the minimized Hamiltonian was only realized where the slope was zero. 




Figure 12. Hamiltonian as a function of u only (From Ross, 2009) 


Pontryagin’s HMC, therefore replaces the Euler-Fagrange equation and 
minimizes H subject to the upper and lower bounds of the control, <u<v^ . 

He also proved that the minimized Hamiltonian is constant with respect to time and that 
for minimum time problems the constant is equal to zero and for fixed time, minimum fuel 
problems is equal to -1. 

It can be shown that from the HMC that a relationship exists between the 
direction of the thrust vector and the direction of the co-state related to velocity. 
Equation (3.10) shows that the angle between them is 180 degrees. 


T 

T 



( 2 . 10 ) 
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D. THE SWITCHING FUNCTION 


Noting that the Hamiltonian as a function of T is linear, and is bounded subject 
toO<r the right hand side of (2.9), which is the coefficient of T in H, becomes 

the switching function, S . 


5;=^cosy0 +—siny0-^ 
m m 


( 2 . 11 ) 


Figure 13 shows the Hamiltonian as a linear function of T, where S is either the 
positive or negative slope of the line. This result implies that when the switching 
function is negative, the thrust is and when positive, the thrust is equal to zero. 

These two cases refer to what is known as bang-bang control, where the control is either 
maximum or minimum. There is a special case where the switching function is equal to 
zero which implies that the thrust lies between the maximum and minimum bounds that is 
covered in the next section. It can be shown that for unconstrained steering controls, the 
HMC yields the same result as the Euler-Lagrange equation, so Pontryagin’s Principle is 
the general case of the Euler-Lagrange condition (Ross, 2009). 




Eigure 13. H as a function of T only (Erom Ross, 2009) 
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E. NECESSARY CONDITIONS FOR SINGULAR ARC OPTIMALITY 


The case where the switching function equals zero occurs when an intermediate, 
or singular, thrust arc appears as part of the optimal trajectory where the thrust will 
assume a value between the upper and lower bounds. Singular optimal control theory 
was developed as a result of the examination of intermediate thrust arcs in central 
inverse-square gravity fields and what resulted was the Higher Order Maximum 
Principle. The Principle follows the results of Pontryagin’s HMC and then takes the 
derivative of the switching function repeatedly to develop the necessary conditions for an 
optimal singular arc. This requires the derivation of an explicit function of the thrust in 
order to develop a singular optimal control law. Since the switching function equals zero, 
each subsequent derivative must clearly also equal zero. Taking the derivative of the 
switching function with respect to time shows that T does not appear explicitly until the 
fourth derivative, which implies that the problem contains a second order singular control 
(Park, C., Yan, H. Gong, Q., and Ross, I.M., 2010). Determining that the fourth order 
derivative is equal to zero is especially important when verifying numerical solutions as 
the results are typically close, but never exactly equal to zero. 

The problem Yan et al. were solving was identical to the problem solved in this 
thesis, the only difference in that the rocket specifications were for the Orion spacecraft 
main engines. Since the equations of motion and dynamics of the problems are identical, 
the derived necessary conditions hold here. 


dS _ 

dt niyjX^-X^ 


( 2 . 12 ) 


d^S 

dt^ 
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dt m^JX^ -X^ dt 


(2.14) 
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(2.15) 


d^S _ 1 d 

myJX^ -X^ dt 

The resulting third and fourth order derivative equations are rather long, and the 
eontrol variable T only appears explieitly at the fourth derivative. The explieit forms for 
the first through fourth derivatives of S with respeet to time are found in Appendix A. 
With these equations, it is now possible to verify the optimality of singular ares that 
appear as part of an optimal earth return trajeetory, however as alluded to previously, 
numerieally validating equality to zero ean be troublesome. Therefore applying another 
result using the ealeulus of variations known as the Generalized Legendre-Clebseh, or 
General Convexity, eondition, see Equation (3.16), the numerieal verifieation is made 
straightforward (Bryson, 1975). In the equation below, the subseript i represents the 
terms for the earth and sun. Terms without a sub-seript are the terms for the moon. 




dT dt^ 


dH ^ 
dT, 


d d^S 
dT dt^ 


> 0 



9Mi[rrK) 15//(r-/l,)' 


l^M (f -^v)' 

ri[K-K) 


(2.16) 


F. TRANSVERSALITY CONDITION 

The derivation of the transversality eonditions are outlined in Appendix B and the 
following eondition holds for the minimum fuel problem. 



(2.17) 
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G. SUMMARY OF NECESSARY CONDITIONS FOR OPTIMALITY 

In summary, the following necessary conditions hold for optimal control: 

1. Hamiltonian is constant. (-Ifor minimum time, 0 for fixed-time, minimum 
fuel) 

2. f = 

T 1, 


3. r = 


0, when 5 > 0 

T^max whenS<0 

Singular when 5=0 


> where S = 

m v„ 


4. Verify optimality of singular arcs. 


a. Numerically verify that —— = 0 for singular arcs. 

dt 


b. Check Generalized Legendre-Clebsch Condition, 


d d^S 
dT dt^ 


> 0 


5. = (for minimum fuel cases only) 


H. GENERIC PROBLEM FORMULATION 

The format for the problems solved in this thesis follow Pontryagin’s method for 
solving optimal control problems. The format is useful, as the formulation of the 
problems in DIDO© are made easier as the boundary conditions and constraints are 
clearly stated. There are primarily two types of optimal controls problems that are 
solved, namely a minimum time and minimum fuel. Each type of problem has a special 
set of boundary conditions that are expressed in the problem formulation. For example, 
the minimum time solution will have a free condition on the final mass and final time. In 
the coding of the problem, the final mass will have to be fixed at zero, or close to zero to 
prevent an unrealistic negative mass solution. In the minimum fuel problem formulation, 
the final time would be fixed and the final mass would be free. 
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Each of the problems solved in this thesis inelude a third physieal dimension, z, 
start with the following format and the eonditions derived for the two-dimensional ease 
deseribed earlier still hold. 


Generie Problem Formulation 
Where v = [v,y,z,v^,v^,v^,m] and u=\T,a,p\, 


Minimize or 

Subjectto X=f{x,u,t^ 

x{t,) = x\x[tf) = x^' 

And any other boundary conditions and/or constraints 
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IV. USING DIDO© 


A. INTRODUCTION 

The use of eomputer-based software paekages sueh as DIDO© allow for solving 
of eomplex optimal eontrol problems, however there are many ehallenges assoeiated with 
solving the types of problems diseussed in this thesis. Primarily, a solid understanding of 
the underlying physies of the problem being solved should be well understood, espeeially 
during the verifieation and validation proeess. This also ineludes an understanding of 
what is being optimized in the problem and eheeking to see if the result makes physieal 
sense as well as meeting the neeessary eonditions for optimality. Finally, an 
understanding of the limits of numerieal eomputation is useful in order to make the eode 
run more effieiently. The verifieation and validation steps taken for the problems in this 
thesis are deseribed, following a brief deseription of how the problems were formulated 
in DIDO©. 

B. PROBLEM FORMULATION IN DIDO© 

DIDO© is a software module that works with Matlab. It implements Legendre 
pseudospeetral methods whieh are used to solve eomplex, non-linear optimal eontrol 
problems (Ross I.M. , 2007). Its primary advantage is that it ean solve the problems with 
as few or many diserete points, or LGL nodes, as defined by the user, whieh makes the 
eomputation time relatively short. However, this requires at iterative proeess of 
inereasing the number of nodes, until the desired fidelity is aehieved. It has a 
straightforward interfaee that allows the user to set up an optimal eontrol problem using 
the previously deseribed Generie Problem Formulation. The problem formulation is 
broken up into two parts, namely the eost funetion, whieh speeifies what design 
parameter is being minimized, and the eonstraints. The eonstraints ean be further divided 
into two parts. The spaeeeraft in eaeh of the problems is eonstrained to move subjeet to 
the equations of motion, or the dynamies. In eaeh of the problems, there is also a set of 
boundary eonditions deseribing the initial and final states of the spaeeeraft whieh were 
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given parameters. The problem formulation in DIDO© follows this separation of the 
problem formulation by requiring separate M-files that are ealled by a main problem file. 
An additional file may be used to speeify a path eonstraint, but it is not always required, 
and is typieally problem dependent. 

The files used in addition to the main program file, then are: 

• Cost Function - outputs are the endpoint eost and running eost. In the 
problems being examined, only endpoint eosts are required and are 
problem speeifie. In the following oases it is either the final time, 
tj- (minimum time ease), or the final mass, (minimum fuel ease). 

• Dynamics Function - provides the differential equation 
x{t) = f[x{t),u{t),t) as desoribed by the problem dynamios. The same 
set of dynamios equations is used for all oases studied. 

• Events Function - used to desoribe the boundary oonditions to inolude the 

initial and final states of the problem. For example, in the minimum-fuel, 
oonstrained-time problem, the initial state is 
Xq = [x°, and the final state is 

Xf = ,y^ ,z^ ,v/ ,v/, vj ] , noting that there is no final mass speoified 

sinoe that is the parameter being solved for. 

The main program file oalls these files and then requires the establishment of 

upper and lower bounds on the states and the oontrols. This is important step for 

primarily oomputational reasons. It provides a range of values that the DIDO© algorithm 

oan use in the output of state and oontrol trajeotories. Given that the physioal 

displaoement and velooities of the spaeeeraft should be limited to a region near the earth 

and moon, the upper bounds of the state variables should be provided aeeordingly. 

However, if the bounds are too small, the resulting trajeetories might result in either an 

infeasible or non-optimal solution. Exeessively large bounds will start to affeet 

eomputation time. DIDO© also allows for an initial guess of the trajeetories that gives 

the program a starting point, however it is possible to get aeeeptable results without an 
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initial guess. The aetual algorithm that solves the problem is ealled in Matlab by a single 
eommand line and its outputs are the final eost, as well as primal and dual struetures. 
The primal strueture ineludes the state and eontrol veetors and assoeiated diserete times 
with their lengths being equal to the number of LGL nodes used. The dual strueture 
ineludes the eo-state and Hamiltonian veetors, also equal to the number of nodes. These 
outputs ean then be used for further analysis for verifieation of optimality and feasibility. 

C. BEST PRACTICES 


1. Scaling 


Most eomputational numerieal methods make use of unit sealing, without loss of 
generality, to inerease eomputational effieieney by inereasing rate of eonvergenee while 
inereasing aeeuraey (Ross I.M., 2009). There are no set ground rules on how to seale 
units, however a uniform range of values ean make the program run smoother (Betts, 
2001). In astrodynamies problems in partieular, sealing beeomes an important faetor in 
problem formulation beeause of the seale at whieh spaeeeraft are operating. Typieally, 
the range of values for the distanees travelled ean be on the order of 10^ meters, but the 
timeframe of interest may only be on the order of a few hours. Also to be eonsidered are 
the values of forees in kilonewtons, for example, whieh depend on the standard units of 
mass, distanee, and time. 

For the problems diseussed in this thesis, the following “eanonieaf’ units were 
used to seale the variables for eomputation. 

• 1 Displaeement Unit (DU) = 1,737 km (mean radius of the Moon) 


• 1 Mass Unit (MU) = 20,339.9 kg (wet mass of spaeeeraft) 


\DU^ 


, .1.033.9sec 


• 1 Time Unit (TU) =-y )| 4,902.8km7see" 

Where is the moon's gravitational eonstant. 


^ \ TT u 1,737km 

• 1 Veloeity Unit =-=-= 1.68 km / sec 

\TU 1,033.9 sec 
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1 Force Unit = 


\MU XDU 
\TU^ 


20,339.9^,.Uyfa. 
(l,003.9sec) 


The units sealed in this fashion are suffieiently elose in order of magnitude to 
maintain eomputational effieieney and are able to be eonverted to and from the aetual 
(un-sealed) values for analysis. 


2. LGL Nodes 


Typieally, these types of eomplex optimal eontrol problems eannot be solved in a 
single step with a large number of LGL nodes primarily due to eomputation time. The 
software used allows the user to ehoose the number of diserete points as the output and 
allows for taking the output of a lower node solution and using it as a guess for a higher 
node solution. If the problems are well sealed and the problem is properly formatted, this 
allows for a rapidly eonverging solution without exeess eomputational time. In the 
problems studied in this thesis, the typieal starting point was at 30 or 40 LGL points and 
then eaeh ineremental step inerease in nodes was of the same order. The exeeption was 
the study of the high thrust singular are in whieh a mueh lower number of nodes was used 
to start beeause of the short time duration of the event. 


3. Providing an Initial Guess 

It is possible, though not required, to provide DIDO© with a “guess” trajeetory 
for the eontrols and states. Providing a reasonable guess will, in most eireumstanees, 
produee a faster result. Eaeh of the eases studied in this thesis started with a low number 
of LGL nodes and with no initial guess, however as the number of nodes and subsequent 
fidelity of the solutions were inereased, the preeeding solutions were used as an initial 
guess for the next iteration. This eoneept of using an initial guess was used in two 
different ways depending on the situation. 

Starting at the low fidelity solution, typieally 30 to 40 LGL nodes, it is possible to 
take the resulting primal strueture (eontrols, diserete time values, and states) and use 
them as the initial guess for a subsequent run of DIDO© at a higher number of nodes. 

Experienee shows that inereasing in inerements of no more that 30 to 40 works best. 
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This method is called bootstrapping, but cannot be used without analysis of the preceding 
result so as not to propagate errors caused by bad problem formulation, scaling, etc. 
Therefore, each step of the problem solution process should be checked for optimality 
conditions, feasibility, and other verification and validation techniques discussed below. 

The difficulty with the angles is that as previously mentioned they are unbounded 
controls, however, in coding the problem it is typically unwise to allow infinite bounds 
for computational reasons. The problem as coded must establish small enough bounds so 
that it can be handled by the computer, but must also be large enough to facilitate the 
possible range of solutions. Witha,y0 e [-2;r,2;r], this allowed for a wide range of 
angles from which the solution could be found. 

The resulting control angles are allowed within the program to reach the limits of 
the bounded region, however when this happens it is difficult to determining whether the 
control has hit a “physical” limit imposed by the code or the angle is equal to that limit, 
which in this case is zero degrees. Because of the non-unique nature of trigonometric 
functions, a one-to-one mapping of solutions is not possible and computer algorithms 
typically cannot discern the difference between them (i.e. sin(;r/2) = sin(5;T/2) = 1). 
Therefore, for any computer algorithm additional conditions may have to be imposed. 

In certain cases, angle reduction formulas could be used to convert the output 
angles into a set of angles that fall between a,/3 [-;r, ;r] which would then provide the 

full range necessary for the problem solution. From Equation (1.16), the Cartesian 
components of the thrust vector can be re-written as: 


=T cos(«)cos(y0) 

=rsin(a)cos(y0) (4.1) 

r=rsm(/?) 


From these equations, explicit formulas for the angles can be derived. 
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1 


(4.2) 


a = tan 


= sin ' 


V 
/ 

V 



r = 




Noting that a singular condition occurs when |F| = 0, this can be negleeted by the 

faet that any eombination of angles has no effeet on a spaeeeraft’s trajeetory when the 
thrust is zero, therefore when this is the ease the angles ean be foreed to zero with no loss 
of generality to the problem. Figures 14 and 15 show how a set of resulting angles from 
DIDO© ean be reealeulated using Equation (4.2) whieh ean then be used as the initial 
guess for the angles for the next iteration. 



Figure 14. Reealeulation of Azimuth Control Angle 
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Figure 15. Recalculation of Elevation Control Angle 

As the problem increases in fidelity and the control angles start to fall within a 
tighter range of values, it is then possible to tighten up the coded bounds to help speed up 
computation. 

D. VERIFICATION AND VALIDATION 
I. Examining Initial Output 

One of the first steps in the analysis conducted was to generate plots of the state 
and control trajectories as well as the evolution of the Hamiltonian and the co-states. One 
of the first checks is to determine if any unexpected or undesired limits of the imposed 
bounds were met. In the case of thrust, the limits are expected because of the actual 
physical limitations of the thruster (no negative thrust and a maximum attainable thrust). 
Hitting the bounds on the control angles was covered above. For a minimum time 
problem, the final time should not hit the upper bound because it usually results in an 
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infeasible solution indieating that the time horizon should be longer. In the minimum 
fuel problem, the upper time bound is typieally fixed, and therefore eoded appropriately. 

2. Checking Necessary Conditions for Optimality 

From the dual variables output of DIDO© the Hamiltonian and eo-state veetors 
ean be plotted with respeet to time and using the neeessary eonditions from Chapter III 
ean be used to verify the optimality of the solution. 

3. Checking Feasibility 

The DIDO© endpoint eonditions ean be validated against the given endpoint 
eonditions, however the purpose of performing the feasibility analysis is to validate the 
DIDO© solution and determine if the solution has enough fidelity to be “flyable.” A 
method of verifying the feasibility of the solution is to take the same initial eonditions for 
the problem and propagate the same dynamies, using Matlab’s ode45 solver, over the 
same timeframe to evaluate how elosely the state trajeetories follow the DIDO© solution. 
This requires an interpolation of the generated eontrols for areas between the nodes sinee 
there many more points used in the ode45 solver eompared to the sparse number of points 
generated by DIDO©. Plotting the DIDO© and propagated state trajeetories 
simultaneously ean provide a good visual indieator of whether the solution is making 
physieal sense, assuming the dynamies are eorreet. However, for problems that eover a 
large time horizon, the results will tend to diverge as time goes on and this is primarily a 
result of interpolation and propagation errors inherent to numerieal solutions. 

A teehnique that has been proven to be sueeessful in redueing the propagation 
errors of the interpolated eontrols from the optimal solution generated by DIDO© is 
ealled the Bellman Pseudospeetral (PS) Method (Ross, I.M., Gong, Q., and Sekhavat, P., 
2008). Its applieation is based on the Bellman Prineiple of Optimality that states that for 
an optimal path between two points, any point taken along the path will follow the 
optimal path to the final point. 

The Bellman Prineiple of Optimality states that along a given optimal trajeetory, a 

segment starting from some intermediate point on that trajeetory and ending at the 
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original termination point will also be optimal. If some trajeetory AB meets the 
necessary conditions for optimality, then a segment starting from point C on the path will 
follow the same optimal path to the final point. In Figure 16 the total cost of the optimal 
trajectory AB is . If the point C lies on the trajectory AB , then the trajectory from CB 

will also be an optimal path on AB , and 7 q = 7j + (Kirk, 1970). 



A 


Figure 16. Bellman Principle 

Yan et al. use the Bellman PS Method as a means for mesh refinement and taking 
advantage of DIDO© node distribution which puts a higher concentration of nodes near 
the beginning and end of the time horizon of the optimal trajectories. Hence, the 
resolution of data near the beginning and end of the trajectory is higher than in areas in 
the middle. Increasing the number of nodes would only offer a limited increase in 
resolution near the middle of the trajectory which would come at a high cost in 
computational efficiency. They found that by applying the technique in this manner, the 
interpolation and propagation errors were significantly mitigated by a few orders of 
magnitude (Yan et ah, 2010). The application of the corollary statement was used in the 
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detailed analysis of the singular are whieh used two points on the trajeetory prior to and 
after the seeond TEI maneuver of the main engine solution as a segment whieh lay on an 
existing optimal path. 

4. Comparison to Previous Studies 

Finally, the results of the ease studies in this thesis ean be eompared to the results 
by Yan et al. study of optimal moon to earth trajeetories using the Orion spaeeeraft main 
engines as well as the study eondueted by Park et al. on numerieal verifieation of singular 
ares for the same problem using a low thrust engine. The roeket engine parameters, 
resulting fuel eonsumption, and eontrol maneuvers ean be eompared for analysis. 


44 



V. AUXILIARY ENGINE EOR MOON TO EARTH RETURN 

MISSION 


For the three eases for the moon to earth return mission, the seenario depieted is a 
loss of the main engine with the same mission time requirement and fuel eonstraint. The 
initial position and veloeity, or state, of the spaeeeraft is at a designated point in a low- 
lunar orbit from whieh the mission eommenees. The final state aehieves the eonditions 
that will take the spaeeeraft to the earth interfaee point without any further bum 
maneuvers. The initial epoeh is given as April 4, 2024, 15:30:00 TDT. Eaeh of the three 
eases studied, the final position and veloeity remain the same and the results of ehanging 
the boundary eonditions on time and fuel eonsumed are examined for the purpose of 
determining the overall feasibility of using the auxiliary engines as an alternative to the 
main engines for the return mission. 

The spaeeeraft eharaeteristies and fixed initial and terminal eonditions are as 
follows (Searritt, S., Marehand, B., and Weeks, M., 2009): 

• Initial Mass: 20,339.9 kg (wet mass) 

• Total Fuel Mass: 8,063.65 kg 

• Auxiliary Engine Thmst: 4,448.0 N 

• Auxiliary Engine Isp: 309 see 

• Initial State in the J2000 Moon eentered inertial frame: 


Xo 


-1,236.7970783385588 km' 

3^0 


1,268.1142350088496 km 

^0 


468.38317094160635 km 

Ao 


0.0329108058365355 km/sec 



0.589269803607714 km/sec 

AzoJ 


-1.528058717568413 km/sec 
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Final State (to reach specified earth interface point) 






2,106.40349727063 km' 

3^/ 


3,596.40440859595 km 



608.167937792834 km 



-0.151018343429892 km/sec 

^yf 


-0.321178999121464 km/sec 



-1.79844521633271 km/sec 


A. THE MINIMUM TIME, FUEL FREE PROBLEM 

I. Problem F ormulation 

The first case of this problem was formulated as a minimum time solution with 
the final mass allowed to be free, but limited to a lower bound of 10% of the total 
spacecraft wet mass. While this problem did not satisfy the fuel requirements, it was an 
important first step in determining an overall feasible engineering solution. By making 
both the time and fuel constraints ‘free,’ it was possible to first determine the feasibility 
of the mission time requirement imposed on this segment of the return mission. With 
time as the cost function to be minimized, the spacecraft was allowed to bum as much 
fuel as needed, even beyond the maximum fuel capacity of 8,063.65 kg. In the Matlab 
code, however, an artificial lower limit of 10% of the spacecraft dry mass was imposed to 
prevent the mass from going negative. Provided an optimal solution existed, the goal was 
to determine if the minimal final time was less than the 48-hour mission constraint. This 
in turn implied that the use of the specified auxiliary engine would be feasible, at least in 
terms of the mission time constraint. 
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Minimum Time Formulation 


Where x = 


X, y, z, v^, v^, u 


, m and u = \T,a,p\, 


Minimize 
Subject to 


< 


i= 

z= 


V, V 

Me{x-x^) 

V,= 3 

3 


fE 

■ ■■ Muy 

ME{y-yE) 

Vy= yu 3 

3 

ru 

^E 

_ “ Mm ^ 

Me{z-Ze) 


3 

fM 

T 

m= - 

^E 

x{t,) = x\x(t^.) = x^ 
tf and are free 

0<T<T 

max 


//j (x - ^ T eos a cos /? 

r/ m 

l^s{y-ys) T sin a cos /? 
r/ m 

^ Tsin^ 
r/ m 


2. Results and Analysis 

This case was solved, starting at 30 LGL nodes and bootstrapping the solution up 
to 80 LGL nodes in order to achieve a smooth and convergent solution. As expected, in 
order to achieve the desired end point condition in minimum time, the spacecraft will 
bum fuel at a maximum rate, using maximum thmst, throughout the maneuver as 
depicted in Figure 17. Figures 18 and 19 show that the azimuth and elevation control 
angles are continuous and while they are not physical constraints on the problem, the 
angles do not hit the bounds specified in the Matlab code. While the fuel mass consumed 
exceeds the actual fuel capacity of the spacecraft, the total time of the maneuver is much 
less than the 48 hours required by the mission. This implies that the auxiliary engines. 


47 



subject to the dynamics of the problem, can be used and will satisfy the time 
requirements of the mission. Further analysis follows to verify that the auxiliary engines 
can be used to meet the fuel requirements as well. 
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Figure 17. Thrust Profde 
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Figure 18. Azimuth Control Profile 



Figure 19. Elevation Control Profile 
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a. Feasibility of the Solution 

As part of the feasibility analysis, the solution is propagated using the 
ode45 funetion in Matlab, using the initial boundary eondition and interpolating the 
controls over the course of the trajectory and comparing the continuous solution against 
the DIDO© solution of discrete points. By examining the profiles of the state vectors for 
position, velocity, and mass, it was found that the DIDO© solution follows the 
propagated solution and is deemed a feasible solution. Figures 20, 21, and 22 show the 
displacement, velocity, and mass trajectories. Table 3 shows the differences between the 
DIDO© and the propagated solutions. While the DIDO© solution achieves the specified 
boundary conditions, the propagated solution shows an error greater than the required 
limits of 1 km for displacement and 1 m/s for velocity. 



Figure 20. Displacement Vector 
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Figure 21. Tangential Veloeity Veetor 



Figure 22. Fuel Consumption 
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Final Position Error (km) 

Final Veloeity Error (m/s) 

DIDO© 

4.5e-13 

2.2e-13 

Propagated Solution 

17.0 

4.2 


Table 3. DIDO vs. Propagated Solution 


The diserepaneies in the final position and veloeities are due to the 
numerieal sensitivities of this type of problem. Very slight interpolation errors at the start 
of the orbit transfer ean eause very large errors at the terminal position (Yan et al., 2010). 
In this partieular problem, the total time horizon is shorter than the 48-hour requirement 
and thus results in a relatively small error. Within two iterations of the Bellman 
teehnique, the required end-point errors are aehieved. For brevity, the Bellman teehnique 
is demonstrated for this problem only, but ean be used to resolve the differenee between 
the DIDO© solution and the propagated solution in the minimum fuel ease. 

The Bellman pseudo-speetral (PS) teehnique is summarized as follows 
(Ross et al., 2008); 

1. Solve the optimal eontrol problem using a reasonably low number 
of LGL nodes to generate a diserete time solution. 

2. Partition the time interval [^tQ,t^^into Bellman segments 

< t' <... < =f. The segments do not need to be uniformly 

spaeed, however extend from a speeified time until the terminal 
time, tf. 

3. Propagate the system dynamies from t°to t'using Tgas the initial 
eondition and any method of eontinuous-time reeonstruetion of the 
eontrols, w, (t),t e J . 

4. Set Tg = v' ^ and = t' and go baek to Step 1. 

This algorithm ends when the final required eonditions are met. 
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Using the Bellman technique in two iterations, the final required accuracy 
was attained. The first Bellman segment [^q,t^Jwas set to start at 0.8363 hours. For 
simplicity, a time associated with one of the discrete time points was chosen. For the 
second Bellman segment J, a start time of 0.9640 hours was chosen. Table 4 shows 

the segment partition and the resulting final position and velocity errors. After two 
iterations, the required errors for both position and velocity are attained. 



Position Error (km) 

Velocity Error (m/s) 

Without Bellman technique 

17.0 

4.2 

First Bellman segment 

1.2 

0.31 

Second Bellman segment 

0.29 

0.16 


Table 4. Effects of the Bellman PS Technique 


Another verification of the feasibility of the solution is to examine the 
osculating orbital elements throughout the maneuver. Figure 23 shows the evolution of 
the changing semi-major axis throughout the maneuver. Towards the end of the 
maneuver, the semi-major axis goes to infinity which is consistent with a parabolic 
trajectory, but then completes the maneuver at a finite semi-major axis value. Figure 24 
shows the eccentricity constantly changing throughout the orbit and eventually passes 
through e = 1, a parabola, and ends on a value greater than one which defines a 
hyperbolic trajectory. Both parabolic and hyperbolic trajectories imply that a spacecraft 
has sufficient energy to overcome the gravity well of the body which it has been orbiting 
(Curtis, 2005). This result shows that the spacecraft in this case has escaped the moon’s 
gravity. Finally, Figure 25 shows a continual change in the inclination of the trajectory 
with respect to the moon’s equator. The final states of the semi-major axis, eccentricity, 
and inclination are consistent with the end point conditions found for the minimum fuel 
solution for the main engines (Yan et ah, 2010). 
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Figure 23. Semi-major Axis 



Figure 24. Eccentricity 
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Figure 25. Inclination 

b. Verification of the Optimality of the Solution 

The necessary conditions for optimality are met and are examined as 
prescribed in Chapter III. In Figure 26, the Hamiltonian is constant throughout the 
maneuver and its value is consistent for a minimum time solution. From the 
transversality condition, the co-state related to mass has a terminal value equal to zero as 
shown in Figure 27. 
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Figure 26. Flamiltonian Evolution 



Figure 27. Co-state (A,m) Evolution 
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The condition derived from the HMC that the angle of thrust is opposite 
the velocity co-state is depicted in Figure 28. During the maneuver, the angle maintains a 
value of 180 degrees as the thrust is “on” throughout. And finally, the Switching 
Function as defined by the HMC should be negative if the constrained control, thrust, is 
at maximum which is true in this case. Figure 29 shows that the switching function is 
indeed negative throughout the maneuver which is consistent with the thrust being at its 
maximum value throughout the maneuver. 



Figure 2 8. Primer Angle 
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Figure 29. Switching Function 

3. Conclusion for Minimum Time Case 

The results show that the solution determined in the minimum time problem is 
optimal with respect to time. The maneuver achieves the final boundary conditions well 
within the 48-hour time constraint for the TEI maneuvers for the moon to earth return 
mission, which informs that the auxiliary engines are capable for meeting the mission 
time requirements. The next step is to determine the minimum fuel required that satisfies 
the mission time requirements and determine the feasibility of using the auxiliary engines 
subject to the fuel capacity constraints of the spacecraft. 


- Thrust (F) 

. Switching Function (S) 
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B. THE MINIMUM FUEL, FIXED TERMINAL TIME PROBLEM 

I. Problem F ormulation 

The auxiliary engine speeifieations as well as the initial and terminal boundary 
eonditions remain the same as in the minimum time optimal eontrol problem formulation. 
In this ease, the problem is formulated as to minimize the fuel eonsumption over a fixed 
time horizon. The time horizon is eonstrained to 48 hours and represents the allowable 
window for the Orion spaeeeraft to conduet the TEI burns for an expedited return to 
earth. The goal was to find a feasible solution to the problem of eondueting the return 
mission using solely the auxiliary engines. Beeause of the roeket eharaeteristies, the 
auxiliary engine uses more fuel for a given duration of burn than does the main engine. 
This solution, if found, was to show that the spaeeeraft would have enough fuel onboard 
to eonduet the return mission. Additionally, verification of the optimality of the solution 
is necessary. 
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Minimum Time Formulation 


Where x = 


X, y, z, v^, v^, u 


, m and u = \T,a,p\, 


Minimize 
Subject to 


< 


i= 

j= 

z= 




'M 


-yu 


MmS f^E{y Se) 


= -Zm - 


I^M 


z ^e{z-Ztz) 


m= - 


T 

v„ 


x{t,) = x\x(tf) = x^ 
tf - to = 48 hrs, m is free 
0<T<T 


//j (x - ^ T eos a cos /? 

r/ m 

l^s{y-ys) T sin a cos /? 
r/ m 

Ms{z-Zs) ^ Tsmp 
r/ m 


2 . Results and Analysis 

Figures 30, 31, and 32 show the resulting control profiles for the minimum fuel 
problem. The problem was initially started using 30 LGL nodes and then bootstrapping 
the solution to a final solution using 300 LGL nodes. The goal was to balance 
computation time by limiting the number of nodes at the same time as achieving a 
smooth and viable control solution. As the number of LGL nodes was increased for each 
iteration of the solution, the control angles were forced to zero whenever the thrust was 
zero. This helped speed up the computation of the DIDO© results and was assumed that 
this would have no effect on the optimal solution since the control angles only influence 
the spacecraft trajectory whenever there is a thrust associated with them. The thrust 
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Thrust (KN 


profile is similar to the main engine trajectory solved by Yan et al. in that it contains three 
distinct TEI maneuvers, with the second bum having the characteristics of a singular arc. 



Figure 30. Thmst Control Stmcture 
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Figure 31. Azimuth Control Angle 



Figure 32. Elevation Control Angle 
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The entire series of maneuvers are completed within the bounded time horizon 
and, as depicted in Figure 33, uses a total of 7,215.6 kg of fuel, less than the given fuel 
capacity of the spacecraft which is 8,063.65 kg. 



Figure 33. Fuel Consumption 


a. Feasibility of the Solution 

A feasibility analysis of the solution was conducted, similarly to the 
previous case, where the solution is propagated using the Matlab ode45 function. The 
controls are interpolated over the propagated dynamics to see how well the discrete 
DIDO© solution maps to the integrated solution. Figures 34 and 35 illustrate that the 
displacement and tangential velocities of the discrete solution follow the propagated state 
trajectories. Some divergence exists towards the end of the trajectory due to numerical 
propagation errors. The difference in the terminal condition of the propagated solution 
and the DIDO© solutions are 502.1 km in displacement and 0.10 km/sec in velocity. The 
variation in final mass was approximately 0.1 kg. While DIDO© finds a solution that 
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meets the given boundary conditions, a resolution to the validation and verification of the 
results need be conducted since the propagated solution does not meet the displacement 
error and velocity error requirements of 1 km and 1 m/s respectively. These errors can be 
eliminated using the Bellman technique described in the minimum time problem. 
However, the initial analysis shows that since the propagated solution and discrete 
solutions are closely matched, the solution appears feasible. 



Figure 34. Displacement State Trajectory 
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Figure 35. Velocity State Trajectory 

Again, like the analysis performed in the minimum time problem, the 
evolution of the osculating orbital elements was investigated. In Figure 36, the semi¬ 
major axis changes primarily due to the first and last TEI bums. The first burn is an 
apoapsis raising maneuver and the final burn provides the kinetic energy required for the 
spacecraft to achieve escape velocity from the moon. In a parabolic trajectory, the semi¬ 
major axis is not defined, hence the singular spike in the figure. Figure 37 shows that the 
eccentricity of the orbit increases to a value greater than one, which means that the 
spacecraft departs the moon’s orbit on a hyperbolic trajectory, requiring a transition 
through a parabolic trajectory where the eccentricity equals exactly one. Figure 38 shows 
an inclination change of about 15 degrees corresponding to the second, singular arc 
maneuver. The terminal values of the osculating orbital elements match those of the 
main engine TEI maneuver sequence (Yan et ah, 2010). It should be noted that in each of 
the TEI bums, a combination of inclination and eccentricity changes is found. This 
appears to be the case since each of the maneuvers is finite in duration, rather than 
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impulsive. This effect is amplified because of the lower thrust engine which requires a 
longer burn time to perform a given maneuver than one with higher thrust. 
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Figure 36. Semi-major Axis With Respect to Moon Centered Frame 


66 






2 


1.8 
1.6 
1.4 
1.2 
tn 1 
0.8 
0.6 
0.4 
0.2 


10 


Hyperbolic Trajectory (e>1) ^ 


A, 


20 30 40 

Time (hrs) 


Figure 37. Eccentricity With Respect to Moon Centered Frame 
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Figure 38. Inclination With Respect to Moon Centered Frame 
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A final check to ascertain the feasibility of the solution is to compare the 
final fuel consumption to the results from several authors who solved the same minimum 
fuel earth return problem using the same initial and final boundary conditions, but used 
different engine characteristics. In Yan et ah, as previously described solve the problem 
using the main engine with a thrust of 33.6 kN and an Isp of 326 seconds. Park et al. 
used a rocket engine with a thrust of 3 kN and an Isp of 326 seconds. The results in 
Table 5 show the resulting fuel consumption of the present problem using only the 
auxiliary engines for the entire return mission compared to each of these cases. The fuel 
consumption value of 7,215.6 kg should be close to the value reached by Park with 
similar engine characteristics and greater than the main engine which with a higher thrust 
can perform the maneuvers with shorter bum durations and greater fuel efficiency 
(greater Isp). 


Engine Configuration 

Isp (sec) 

LGL Nodes 

Fuel Consumed (kg) 

Auxiliary Engines (4.4 kN) 

309 

300 

7,215.6 

Auxiliary Engines (3 kN) 

326 

160 

7,245.1 

Main Engines (33.6 kN) 

326 

500 

6,657.3 


Table 5. Fuel Consumption Comparison 


b. Verification of the Optimality of the Solution 

Figures 39 and 40 show the Hamiltonian constant at a value of zero, and 
the transversality condition of the co-state related to mass equal to -1 at the terminal 
point. 
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Figure 3 9. Flamiltonian Evo lution 



Figure 40. Co-state (A,m) Evolution 
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Figure 41 shows the evolution of the switching function over the entire 
maneuver sequence. The switching function is negative for maximum thrust and positive 
for zero thrust, however at the singular arc the switching function also appears to be 
negative. The singular arc is demonstrated to be optimal by examining the derivatives of 
the switching function up to the fourth derivative. Figure 42 shows that the fourth 
derivative is equal to zero during the singular arc. The generalized Legendre-Clebesch 
convexity condition during the singular arc burn is shown in Figure 43. Figure 44 is 
included solely to demonstrate that it is possible to have regions where the Legendre- 
Clebesch condition is non-positive. 
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Figure 41. Thrust and Switching Function 


70 

















Generalized Legendre-Clebesch (scaled) t & d^s/dt^ (scaled) 



Time (hrs) Time (hrs) 


ure 42. 


Derivatives of the Switching Function During Singular Arc 



Figure 43. Generalized Legendre-Clebesch Condition During Singular Arc 
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Figure 44. Generalized Legendre-Clebesch Condition Over Entire Trajectory 
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Figures 45, 46, and 47 show in detail the features of the TEI bums; they 
show the thmst and control angles and also that for each burn, the primer angle between 
the thmst and the velocity co-vector is equal to 180 degrees as necessary to meet the 
Hamiltonian Minimization Condition. With the exception of the sparse nature of the 
control trajectory of the second, singular TEI bum, the control solution appears smooth 
enough for a control solution. From a practical control perspective, however, it is desired 
to increase the fidelity of the discrete solutions into a more continuous solution. 
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Figure 45. First TEI Burn 



Figure 46. Second TEI Bum 
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Figure 47. Third TEI Bum 
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As pointed out in the feasibility analysis, the final displacement and 
velocity errors are due to the interpolation errors propagated over a long time horizon. 
Figures 48 and 49 show how the error between the DIDO© solution and the interpolated 
solution increases over the entire solution trajectory. Using the Pseudospectral Bellman 
Technique used by Ross et al., and demonstrated in the minimum time problem, it is 
possible to minimize the errors caused by interpolation of the controls. 


74 






5 V (km/sec) 



Figure 48. Displacement Error 



Figure 49. Velocity Error 
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c. 


CONCLUSIONS FOR USE OF AUXILIARY ENGINES FOR RETURN 
MISSION 


The resulting optimal eontrol trajeetory shows that in the event of a total main 
engine failure, the auxiliary engines eould be used for the entire return mission, meeting 
the 48-hour minimum time requirement to eomplete the third TEI maneuver. 
Additionally, the fuel optimal solution shows that 7,215.6 kg is the theoretieal minimum 
required fuel to be earried onboard, not ineluding any design margins. This result then 
potentially reduees the total amount of fuel required, thus freeing up additional mass for 
other aspeets of the spaeeeraft design. 


76 



VI. MAIN ENGINE SINGULAR ARC ANALYSIS 


The Yan et al. results for the minimum fuel solution for the earth return mission 
using the main engines, they found a singular arc as part of the control trajectory for the 
second TEI maneuver. In their solution, the singular bum is made up of only three 
distinct, non-zero thmst points. Since a high fidelity control solution is desired, it would 
be advantageous to have additional points from which to approximate a continuous 
solution. As a secondary investigation, NASA is interested in determining the feasibility 
of using the auxiliary engines in place of the main engines for the singular arc TEI 
maneuver. Therefore, the purpose of the singular arc study is to: 


(1) Generate a higher fidelity control solution from the existing singular arc 
data generated by Yan et al. 

(2) Determine the feasibility of using the auxiliary engines for the singular arc 
and ascertain the fuel penalty if any 

(3) Verify and validate the necessary conditions for optimality for any feasible 
solution(s) generated 


The separation between the nodes from the original main engine solution is due 
primarily to its location along the entire trajectory solution. The output data points from 
DIDO© are non-uniformly spaced, being more dense at the beginning and end of the 
trajectory and sparse in the middle. The singular arc was in the sparse region of the 
original solution covering a time span of approximately 29 minutes. The data points 
around the singular arc were spaced approximately 7 minutes apart as compared to the 
points around the final TEI bum had a spacing range of approximately 0.4 to 0.8 minutes 
between them. The approach therefore was to examine the singular arc by fixing 
boundary conditions close to and around the start and completion of the maneuver. Since 
the goal is to increase the number of discrete points for a solution with higher resolution. 
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the problem uses a small number of LGL nodes to eompute a solution and then using the 
bootstrap teehnique to inerease the nodes until the eontrol trajeetories eonverge in a 
smooth fashion. 

Reealling that the Bellman Prineiple of Optimality states that along a given 
optimal trajeetory, a segment starting from some intermediate point on that trajeetory and 
ending at the original termination point will also be optimal. It says nothing about an 
intermediate segment. Suppose there are two points, C and D, forming a segment along 
the same optimal trajeetory. Even though the original optimal trajeetory ean pass through 
these points, it does not mean that the original path on AB between them is optimal on 
the path between C and D . There exists an optimal path CD that may not lie on AB as 
shown in Figure 50. If the minimized eost along AB is , and the sum of the segments 

on AB is J Q = J y + J 2 + J then the eost of an optimal trajeetory CD not on AB must be 

J 2 > J 2 ■ Put another way, 7j + Tj + 3 ^ 1 + 2 + -^3 • 



A 


Figure 50. Bellman Curve and Segment 
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The singular arc generated from the Yan et al. solution can be said to lie on a 
segment of a fuel-optimal trajectory. In order to maintain the same cost along that 

segment, so that from the example above, , then the dashed line trajectory must lie 

on the solid line CD. This implies that the segment has fixed state and time horizon 
boundaries and the control and state trajectories will lie on the original trajectory. In the 
problem formulation in DIDO© this poses a difficulty since to solve an optimal control 
problem, one needs to determine the cost to be minimized. If both mass (part of the state) 
and time are fixed boundary conditions, then there is no cost to minimize. Therefore, the 
problem was formulated to minimize the final fuel cost, fixing the initial boundary 
condition at some Xq and the final boundary condition at 

some Xf =[x'',y^,z^,v/,Vy^,v/] such that the initial and final states are relatively 

close to the initiation and termination of the maneuver. This problem is solved with the 
knowledge that an optimal solution can be found on the segment, but that the cost can be 
no less than the segment on the original optimal path. 

A total of six experiments were run, using three different sets of boundary 
conditions each associated with a different time horizon around the original singular arc 
(THl, TH2, and TH3). For each specified time horizon, DIDO© was supplied with two 
initial guesses to choose from; (1) the original, coarse solution bound by the time 
horizon, and (2) a two-point guess using only the initial and final nodes of the original 
primal solution bounded by the time horizon. Figure 51 shows the thrust profile for the 
singular arc in the original main engine solution and the three time horizons examined. 
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Figure 51. Three Time Horizons Around Original Main Engine Singular Arc 


The spacecraft main engine characteristics (Scarritt et ah, 2009) and fixed initial 
and terminal conditions are as follows: 

• Main Engine Thrust: 33,631.6621 N 

• Main Engine Isp: 326 sec 


• Initial and terminal boundary conditions for THl 
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• Initial and terminal boundary conditions for TH2 


Vo 


-3,098.602243204757 km^ 


^f 


-2,995.8032622945048 km' 

3^0 
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• Initial and terminal boundary conditions for TH3 


Vo 


' -3,066.151333338768 km' 
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Theoretically, the initial mass for each of time horizon should be the same, since 
from the original solution the thruster is not firing, thus not burning fuel at the first three 
nodes of TH3. The differences come from the sensitivity of the scaled mass values. The 
scaled initial mass value for TH3 starts to differ at the fourth significant digit and results 
in a difference of about 1.6 kg. 

• THl and TH2 Initial Mass = 0.825953636218675 

• TH3 Initial Mass = 0.825875239089155 
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A. 


PROBLEM FORMULATION 


The problem formulations for each time horizon are identical, with only the initial 
and terminal boundary conditions changed for each case. For each case, two different 
initial guesses were provided, which in each case influenced the final solution. 

Minimum Time Formulation 
Where x = ^x, y, z, v^, v^, v^, m] and u = [T,a,/3], 
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B. THl RESULTS 

Solving this problem using the THl boundary conditions and supplying two 
different initial guesses to DIDO© results in two different solutions, each satisfying the 
necessary conditions for optimality. The solutions are each solved up to 80 LGL points 
each. Comparing the two cases, the control and state trajectories differ between the fixed 
boundary conditions, however both use the same amount of fuel. Figures 52, 53, and 54 
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show the resulting control trajectories compared to the initial guess trajectory provided. 
Note that the coarse guess is the same trajectory given by the original problem using the 
main engine. Also, the differences in the fuel consumed are negligible as shown in 
Table 6. The solutions for each of the time horizons are compared against the original 
coarse solution for the corresponding time horizon. The fuel consumed for each region 
has differences due to numerical variations of the scaled mass in the original coarse 
solution, therefore the fuel consumed in each of the regions appear to have different 
values. In scaled units the differences are small between the time horizons. 



Fuel Consumed (kg) 

Variation 

Original Coarse Solution 

1,360.2 

— 

Using Coarse Guess 

1,361.0 

0.0006 % 

Using 2-point Guess 

1,361.2 

0.0007 % 


Table 6. THl Fuel Consumption 



Figure 52. Thrust Trajectories for THl 
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Figure 53. Azimuth Angle Trajectories for THl 



Figure 54. Elevation Angle Trajectories for THl 
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Of note is the aliasing that oeeurs when the initial coarse solution is used as an 
initial guess (Ross, I.M., Gong, Q., and Sehkavat, P., 2007). The higher fidelity solution 
follows closely with the original solution. Similarly, the state trajectories are also closely 
aliased when given the original solution. Using the two-point guess, the state and control 
trajectories deviate from the original solution; however, the final minimized cost is the 
same. Therefore, both of these solutions appear feasible and as described below they 
both meet the necessary conditions for optimality. 

Figures 55, 56, and 57 show that the Hamiltonian is constant, the end point 
transversality condition is met, and the primer angle is equal to 180 degrees and thus, the 
Hamiltonian Minimization Conditions are met. Additionally, both solutions meet the 
necessary conditions for optimality for singular arcs and the Legendre-Clebesch 
condition as shown in Figures 58-63. 



Figure 5 5. Hamiltonian Evolution (TH1) 


85 




















0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 

Time (hrs) 












Figure 62. Generalized Legendre-Clebeseh Condition for Coarse Guess Solution 



Figure 63. Generalized Legendre-Clebeseh Condition for 2-Point Guess Solution 
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Examination of the state veetor trajeetories shows that eaeh of the two cases for 
the THl bounds start and complete the maneuver at the specified state boundary 
conditions. While the final mass was not posed as a constraint on the problem, the 
resulting final mass is the same as in the original main engine problem within very tight 
tolerances. In Figures 64 and 65, the position and velocity state vector trajectories for the 
THl singular arc problem are shown and are described below. The position, or 
displacement, state trajectories differ slightly, however the velocity state trajectories are 
more varied due to the immediate relationship between thrust and change in velocity. 




13.9 


-•-Original 



^ 13.85 


Coarse Guess 
2-pt Guess 


13.45 


38.1 38.15 38.2 38.25 38.3 38.35 38.4 38.45 38.5 
Time (hrs) 


Figure 64. Displacement State (r) Trajectories for THl 


90 











Figure 65. Velocity State (v ) Trajectories For THl 


The overall change in mass is the same for the two feasible solutions as shown in 
Figure 66. The mass state trajectory is clearly related to the shorter duration, but higher 
magnitude of thrust in the 2-point guess solution as compared to the solution from the 
coarse guess. 
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Figure 66. Fuel Consumption for THl 


C. TH2 RESULTS 

The same analysis conducted for the TFll solutions was performed for TH2 and 
TH3. The Hamiltonian Minimization Conditions for optimality were met, as well as the 
additional conditions for singular arcs. Each case resulted in different feasible solutions 
for the controls, with the state variable trajectories satisfying the fixed boundary 
conditions. And while the state trajectories vary slightly between the cases, the fuel 
consumption is the same between the original coarse solution and the optimal solutions. 
These results further indicate the range of possibilities for conducting the singular arc 
maneuver. 

In THl, the maximum thrust exceeded the maximum thrust of the auxiliary 
engines; however the resulting control trajectories for TH2 and TH3 achieved maximum 
thrust less than the auxiliary engines (4.4 kN). This gives some indication that the 


92 









auxiliary engines can indeed be used for the singular burn maneuver; however the fuel 
consumption will necessarily be greater due to the lower rocket specific impulse as 
compared to the main engines. 

Figures 67, 68, and 69 show the control trajectories for TH2. The maximum 
thrust for both solutions is about the same, around 4 kN. The control angles during the 
bum are identical to the angles arrived at in the THl solutions. They are also found to be 
equal in the TH3 solutions as well. Clearly, a rocket engine with maximum thmst of 4.4 
kN would be able to execute these maneuvers. The difference will be in fuel 
consumption alone. 



Figure 67. Thmst Trajectories for TH2 
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Figure 68. Azimuth Angle Trajectories for TH2 



Figure 69. Elevation Angle Trajectories for TH2 
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Displacement State Trajectory (scaled) 


Neither of the state trajectories have the same aliasing that was seen in THl, 
however the two solutions, using different initial starting guesses, are very similar as seen 
by the control trajectories above and the state trajectories shown in Figures 70, 71 and 72. 



Figure 70. Displacement State (r) Trajectories for TH2 
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Figure 71. Velocity State (v ) Trajectories for TF12 



Figure 72. 


Fuel Consumption for TH2 
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The fuel eonsumption for TH2 is listed in Table 7 and the differenees are again 
negligible as eompared to the original eoarse solution within the speeifie time horizon. 



Fuel Consumed (kg) 

Variation 

Original Coarse Solution 

1,339.7 

— 

Using Coarse Guess 

1,339.1 

0.0007 % 

Using 2-pomt Guess 

1,339.2 

0.0004 % 


Table 7. TH2 Fuel Consumption 


D. TH3 RESULTS 

By expanding the time horizon to TH3, the eontrol trajeetories begin to deviate 
from the original solution. These solutions also meet the optimality eonditions for 
singular ares and are also feasible solutions. In both solutions arrived at from different 
initial guesses, the maximum thrust is again less than 4.4 kN indieating that this partieular 
eontrol profile could be used with the auxiliary engines. The solution from the 2-point 
guess is interesting in that the thrust profile is starting to appear as a bang-bang 
maneuver, however in this situation remains a singular arc solution. In general, these 
solutions are not as “clean” from a control perspective because neither solution has a zero 
to zero thrust trajectory. In both cases there is a non-zero initial thrust as seen in Figure 
73. In Figures 74 and 75 the control angles are approximately the same for both initial 
guesses even though the coarse gain has a single bum while the 2-point guess solution 
has a two burn solution. 
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Figure 73. Thrust Trajectories for TH3 



Figure 74. Azimuth Angle Trajectories for TH3 
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Figure 75. Elevation Angle Trajeetories for TH3 


The state trajeetories as shown in Figures 76, 77, and 78 illustrate a more extreme 
example of varying state trajectories with the minimized cost of fuel being equal, 
emphasizing the result that multiple local minimal solutions exist around the singular arc. 
Table 8 shows the negligible difference in fuel from the original coarse solution for the 
given time horizon. 
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Displacement State Trajectory (scaled) 



Figure 76. Displacement State Trajectories (TH3) 



Figure 77. Velocity Trajectories (TH3) 
100 






















Figure 78. Fuel Consumption (TFI3) 



Fuel Consumed (kg) 

Variation 

Original Coarse Solution 

1,345.8 

— 

Using Coarse Guess 

1,346.3 

0.0007 % 

Using 2-point Guess 

1,346.7 

0.0004 % 


Table 8. Fuel Consumption for TF13 


E. FEASIBILITY OF USING AUXILIARY ENGINES FOR SINGULAR ARC 
MANEUVER 

In the analysis of the singular arc using the main engine, it was found that for the 
TH2 and TFI3 time horizons, it appeared feasible to use the auxiliary engines in place of 
the main engines to conduct the singular arc maneuver. However, the THl solutions 
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using the main engine resulted in thrust profiles that exceeded the capacity of the 
auxiliary engines. It was expected that the lower Isp would result in a higher fuel 
consumption as compared to the main engine, given the same magnitude of thrust. The 
same six simulations as examined in the previous section were performed using the same 
optimal control problem formulation, but using the thrust and Isp parameters for the 
auxiliary engines. 

As expected, the control profiles around TH2 and TH3 using the auxiliary engines 
were identical to those found using the main engines. The THl solution for the auxiliary 
engines however, was constrained in maximum thrust and resulted in a different, yet 
feasible control structure. Figure 79 shows the different thrust profiles attained in the 
THl time horizon using the same initial primal guess from the main engine singular arc. 
In this case, the two different initial guesses resulted in very similar control profiles. 



Figure 79. Auxiliary Engine Thrust Profile (THl) 
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In all three time horizons, using the auxiliary engines resulted in an increase in 
fuel consumption as shown in Table 9. The results show that for a relatively small 
increase of fuel, the singular arc maneuver can be conducted using the auxiliary engines. 
In order to solve for the total trajectory fuel consumption where the main engines are 
used for the first and third TEI burn and the auxiliary engine is used for the singular arc 
maneuver, another step using the minimum fuel problem formulation needs to be 
conducted between the boundary conditions of the singular arc time horizon and the 
original boundary conditions provided to solve the full moon-to-earth trajectory problem. 



Main Engine 

(33.1 kN) 

Aux Engine 

(4.4 kN) 

Singular Arc on Full Trajectory (Yan et ah, 2010) 

1,333.6 kg 

— 

THl 

1,361.0 kg 

1,432.6 kg 

TH2 

1,339.1 kg 

1,409.6 kg 

TH3 

1,346.3 kg 

1,417.1 kg 


Table 9. Fuel Consumption Comparison Between Main Engine and Auxiliary Engines 

Over Singular Arc 


F. CONCLUSIONS 

The reason that multiple feasible and optimal solutions appear is because of the 
existence of the singular arc. Recalling that the switching function equals zero for a 
singular arc, it is clear that the resulting magnitude of thrust can lie between the minimum 
and maximum thrust. If the initial and final boundary conditions are established at 
varying positions in time before and after the singular arc, it is possible to generate 
different control solutions and thus, different state trajectories that are also optimal 
solutions. This gives rise to the possibility of multiple feasible solutions between two 
given boundary conditions with a singular arc between them. From an engineering point 
of view, this offers greater flexibility in designing a solution that closer resembles the 
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imposed physical constraints on the available controllers. The difficulty though, is in 
finding a fuel optimal solution that meets the given constraints and requirements from a 
possible myriad of solutions. For example, while some modern rocket engines may have 
the ability to throttle, many typical rockets are forced to apply a bang-bang type control. 
In seeking for a bang-bang solution to perform the same maneuver as the singular arc, 
additional constraints will have to be imposed, thus changing the problem formulation. It 
is possible to use pseudo-spectral tools such as DIDO© to find a range of optimal bang- 
bang solutions, however it is likely that the fuel cost will be greater. 

The singular arc maneuver results in a slight change in the orbits eccentricity 
(semi-major axis); however it appears that the maneuver is used primarily to change the 
orbit plane (inclination). Further study of orbit maneuvers may show that singular arcs 
will appear in fuel optimal orbit transfers whenever an inclination change occurs, 
assuming that the rocket thrust is allowed the full range from null to maximum. To 
further examine the nature of fuel-optimal plane change maneuvers, it would be prudent 
to examine pure plane change maneuvers as well as combined plane change and apogee 
raising (or lowering) maneuvers using similar n-body dynamics equations as used here. 
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APPENDIX A. SINGULAR ARC NECESSARY CONDITIONS 


From Park et al, 2010: 
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Fourth derivative 
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Legendre-Clebesch Condition (convexity) 



d d^S 
dT dt" 


> 0 




APPENDIX B. END-POINT TRANSVERSALITY CONDITION 


The transversality condition states; 

The variable, v, represents the Lagrange multiplier vector associated with the end¬ 
point state condition,The end-point cost . The given end-point 

boundary condition is defined as , vj , v/, v/ ] . As can be seen in each 

of the problems, the only additional information yielded is 

Minimum Fuel Problem End-point Transversality Condition 
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Minimum Time Problem End-point Transversality Condition 


E[v,x[t^)) = tf+v\x^-xf) + v,[y^ 

n(v^-v/) + v^5K- 


yiy (tf) 

) 

K{^f) 


dE 

8Xf 

dyf 

dE 


= v, 


= v^ 


dz 


■ = v^ 


f 

dE 


dv 


■ = Va 


-V/ 

dE 


dE 


■ = v^ 


dv 


■ = Va. 


dE 


dm 


7 


-y^)+Vi[zf-z^) + 

V) + ^6K-v/) 


no 



LIST OF REFERENCES 


Betts, J. T. (2001). Practical methods for optimal control using nonlinear programming. 
Philadelphia; Soeiety for Industrial and Applied Mathematies. 

Bryson, A. E. (1975). Applied optimal control: Optimization, estimation, and Control. 
New York: Hemisphere Publishing Corporation. 

Curtis, H. D. (2005). Orbital mechanics for engineering students. Oxford; Elsevier Etd. 

Jet Propulsion Eaboratory. (n.d.). Horizons system. Retrieved Oetober 29, 2010, from Jet 
Propulsion Eaboratory; http;//ssd.jpl.nasa.gov/?horizons 

Kirk, D. E. (1998). Optimal control theory: An introduction. Mineola; Dover Publishing, 
Ine. 

Park, C., Yan, H. Gong, Q., & Ross, I. M. (2010, Eebruary). Numerieal verifieation of 
singular ares on optimal moon-to-earth transfer. Doeument AAS 10-104. Paper 
presented at 20th AAS/AIAA Space Flight Mechanics Meeting, San Diego. 

Park, C., Yan, H. Gong, Q., & Ross, I. M. (2010, Eebruary). Optimality of singular ares 
for a roeket trajeetory of mayor’s variational problem. Doeument AAS 10-106. 
Paper presented at 20'’^ AAS/AIAA Space Flight Mechanics Meeting, San Diego. 

Pineh, E. R. (1993). Optimal control and the calculus of variations. New York: Oxford 
University Press Ine. 

Ross, I. M. (2007). “A beginner’s guide to DIDO; A MATEAB applieation paekage for 
solving optimal eontrol problems.” Elissar Technical Report TR-711 . Monterey, 
CA. 

Ross, I. M. (2009). A primer on pontryagin’sprinciple in optimal control. San Eraneiseo; 
Collegiate Publishers. 

Ross, I. M., Gong, Q., & Sehkavat, P. (2007). Eow-thrust high-aeeuraey trajeetory 

optimization. Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4, pp. 
921-933. 

Ross, I. M., Gong, Q., & Sehkavat, P. (2008, August). The bellman pseudospeetral 
method. Doeument AIAA 2008-6448. Paper presented SiiAIAAJAAS 
Astrodynamics Specialist Conference and Exhibit, Honolulu. 


Ill 



Scarritt, S., Marchand, B., & Weeks, M. (2009, August). An autonomous onboard 
targeting algorithm using finite thrust maneuvers. Paper presented at AlAA 
Guidance, Navigation, and Control Conference and Exhibit, Chieago. 

Vallado, D. A. (1997). Fundamentals of astradynamics and applications. New York: 
McGraw-Hill. 

Yan, H., Gong, Q., Park, C., Ross, I. M., & Souza, C. (2010, August). High-accuracy 

moon to earth escape trajectory optimization. Paper presented at A/AA Guidance, 
Navigation, and Control Conference, Toronto. 


112 



INITIAL DISTRIBUTION LIST 


1. Defense Teehnieal Information Center 
Ft. Belvoir, VA 

2. Dudley Knox Library 
Naval Postgraduate Sehool 
Monterey, CA 

3. Head, Information Operations and Spaee Integration Braneh 
PLI/PP&O/HQMC 

Washington, DC 

4. Dr. HuiYan 

Dept, of Applied Math and Stats 
MS SOE2 

University of California, Santa Cruz 
Santa Cruz, CA 

5. Dr. Qi Gong 

Dept, of Applied Math and Stats 
MS SOE2 

University of California, Santa Cruz 
Santa Cruz, CA 

6. Omar Wheatley 
Author 
Denver, CO 


113 



